{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model Predictive Control: Aircraft Model\n",
    "\n",
    "RMM, 13 Feb 2021\n",
    "\n",
    "This example replicates the [MPT3 regulation problem example](https://www.mpt3.org/UI/RegulationProblem)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import control as ct\n",
    "import numpy as np\n",
    "import control.optimal as obc\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# model of an aircraft discretized with 0.2s sampling time\n",
    "# Source: https://www.mpt3.org/UI/RegulationProblem\n",
    "A = [[0.99, 0.01, 0.18, -0.09,   0],\n",
    "     [   0, 0.94,    0,  0.29,   0],\n",
    "     [   0, 0.14, 0.81,  -0.9,   0],\n",
    "     [   0, -0.2,    0,  0.95,   0],\n",
    "     [   0, 0.09,    0,     0, 0.9]]\n",
    "B = [[ 0.01, -0.02],\n",
    "     [-0.14,     0],\n",
    "     [ 0.05,  -0.2],\n",
    "     [ 0.02,     0],\n",
    "     [-0.01, 0]]\n",
    "C = [[0, 1, 0, 0, -1],\n",
    "     [0, 0, 1, 0,  0],\n",
    "     [0, 0, 0, 1,  0],\n",
    "     [1, 0, 0, 0,  0]]\n",
    "model = ct.ss(A, B, C, 0, 0.2)\n",
    "\n",
    "# For the simulation we need the full state output\n",
    "sys = ct.ss(A, B, np.eye(5), 0, 0.2)\n",
    "\n",
    "# compute the steady state values for a particular value of the input\n",
    "ud = np.array([0.8, -0.3])\n",
    "xd = np.linalg.inv(np.eye(5) - A) @ B @ ud\n",
    "yd = C @ xd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# computed values will be used as references for the desired\n",
    "# steady state which can be added using \"reference\" filter\n",
    "# model.u.with('reference');\n",
    "# model.u.reference = us;\n",
    "# model.y.with('reference');\n",
    "# model.y.reference = ys;\n",
    "\n",
    "# provide constraints on the system signals\n",
    "constraints = [obc.input_range_constraint(sys, [-5, -6], [5, 6])]\n",
    "\n",
    "# provide penalties on the system signals\n",
    "Q = model.C.transpose() @ np.diag([10, 10, 10, 10]) @ model.C\n",
    "R = np.diag([3, 2])\n",
    "cost = obc.quadratic_cost(model, Q, R, x0=xd, u0=ud)\n",
    "\n",
    "# online MPC controller object is constructed with a horizon 6\n",
    "ctrl = obc.create_mpc_iosystem(model, np.arange(0, 6) * 0.2, cost, constraints)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<InterconnectedSystem>: sys[5]\n",
      "Inputs (2): ['u[0]', 'u[1]']\n",
      "Outputs (5): ['y[0]', 'y[1]', 'y[2]', 'y[3]', 'y[4]']\n",
      "States (17): ['sys[3]_x[0]', 'sys[3]_x[1]', 'sys[3]_x[2]', 'sys[3]_x[3]', 'sys[3]_x[4]', 'sys[4]_x[0]', 'sys[4]_x[1]', 'sys[4]_x[2]', 'sys[4]_x[3]', 'sys[4]_x[4]', 'sys[4]_x[5]', 'sys[4]_x[6]', 'sys[4]_x[7]', 'sys[4]_x[8]', 'sys[4]_x[9]', 'sys[4]_x[10]', 'sys[4]_x[11]']\n",
      "\n",
      "Update: <function InterconnectedSystem.__init__.<locals>.updfcn at 0x167dff0a0>\n",
      "Output: <function InterconnectedSystem.__init__.<locals>.outfcn at 0x167dff130>\n"
     ]
    }
   ],
   "source": [
    "# Define an I/O system implementing model predictive control\n",
    "loop = ct.feedback(sys, ctrl, 1)\n",
    "print(loop)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Computation time = 28.414 seconds\n"
     ]
    }
   ],
   "source": [
    "import time\n",
    "\n",
    "# loop = ClosedLoop(ctrl, model);\n",
    "# x0 = [0, 0, 0, 0, 0]\n",
    "Nsim = 60\n",
    "\n",
    "start = time.time()\n",
    "tout, xout = ct.input_output_response(loop, np.arange(0, Nsim) * 0.2, 0, 0)\n",
    "end = time.time()\n",
    "print(\"Computation time = %g seconds\" % (end-start))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.66523705,  0.01149905,  0.23159795,  0.03076594,  0.00674534])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnYAAAHWCAYAAAD6oMSKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAA9hAAAPYQGoP6dpAABwHElEQVR4nO3dd3wUdf4G8Ge2lySbsiGFNIIUKSqCBUQBPT3FeoqejaKIIiByyp0F/YF6iGI96VgARQU824m9YFdAigpIk1TS2ya7m60zvz8mWRJIFgKbTHbzvH3Na3ZnPzv7yRp2n0z5jiBJkgQiIiIiCnsqpRsgIiIiotBgsCMiIiKKEAx2RERERBGCwY6IiIgoQjDYEREREUUIBjsiIiKiCMFgR0RERBQhGOyIiIiIIgSDHREREVGEYLAjImpi8eLFWLlyZYe8ltPpxJw5c/D11193yOsRUeRjsCMiaqKjg90jjzzCYEdEIcNgR0RERBQhGOyIKOx9//33uOCCCxAdHQ2TyYRhw4bhww8/DDw+Z84cCIJwxPNWrlwJQRCQm5sLAMjKysLOnTvxzTffQBAECIKArKwsAMDXX38NQRCwevVq3HPPPUhOTobRaMSIESOwbdu2ZusdOXIkRo4cecTrTZgwIbC+3NxcJCYmAgAeeeSRwOtNmDABAFBeXo7bb78d6enp0Ov1SExMxDnnnIMvvvjixN4sIopoGqUbICI6Ed988w0uvPBCnHLKKXj55Zeh1+uxePFiXH755XjzzTfx97///ZjX9e6772LMmDGwWCxYvHgxAECv1zerefDBB3H66afjpZdegs1mw5w5czBy5Ehs27YN2dnZx/xaKSkp+OSTT3DxxRdj4sSJuO222wAgEPbGjh2LrVu3Yu7cuejduzdqamqwdetWVFZWHvNrEFHXw2BHRGHt/vvvR1xcHL7++mtERUUBAC677DKcdtppmDlzJq677rpjXtegQYNgNBoRExODs88+u8WaxMREvPvuu4EtgMOHD0evXr0wb948vPjii8f8Wnq9HoMHDwYApKWlHfF6P/zwA2677TZMmjQpsOzKK6885vUTUdfEXbFEFLYcDgc2btyIMWPGBEIdAKjVaowdOxaFhYXYs2dPSF/zxhtvbLZbNzMzE8OGDcOGDRtC+jpnnnkmVq5ciX//+9/4+eef4fV6Q7p+IopMDHZEFLaqq6shSRJSUlKOeCw1NRUAQr7rMjk5ucVloX6dtWvXYvz48XjppZcwdOhQxMfHY9y4cSgpKQnp6xBRZGGwI6KwFRcXB5VKheLi4iMeKyoqAgBYrVYYDAYAgNvtblZTUVHR5tdsKViVlJQgISEhcN9gMBzxWm19PavViueffx65ubnIy8vDvHnz8M477wROriAiagmDHRGFLbPZjLPOOgvvvPMO6uvrA8tFUcTq1auRlpaG3r17B85E/e2335o9/4MPPjhinXq9vtm6Dvfmm29CkqTA/by8PPz444/NzoLNysrC3r17m4W7yspK/Pjjj0e8FoCgrwcAGRkZmDZtGi688EJs3bo1aC0RdW08eYKIwtq8efNw4YUXYtSoUZg5cyZ0Oh0WL16MHTt24M0334QgCBg9ejTi4+MxceJEPProo9BoNFi5ciUKCgqOWN/AgQOxZs0arF27FtnZ2TAYDBg4cGDg8bKyMvztb3/DpEmTYLPZMHv2bBgMBjzwwAOBmrFjx2LZsmW4+eabMWnSJFRWVmL+/PmIiYlp9lrR0dHIzMzE+++/jwsuuADx8fGwWq2Ii4vDqFGjcOONN6Jv376Ijo7G5s2b8cknn+Dqq69uvzeTiMKfREQU5r777jvp/PPPl8xms2Q0GqWzzz5b+uCDD5rVbNq0SRo2bJhkNpul7t27S7Nnz5ZeeuklCYCUk5MTqMvNzZUuuugiKTo6WgIgZWZmSpIkSRs2bJAASK+99po0ffp0KTExUdLr9dK5554r/fLLL0f0tGrVKunkk0+WDAaD1K9fP2nt2rXS+PHjA+tr9MUXX0iDBg2S9Hq9BEAaP3685HK5pMmTJ0unnHKKFBMTIxmNRqlPnz7S7NmzJYfDEeq3j4giiCBJTfYpEBFRi77++muMGjUKb731FsaMGaN0O0RELeIxdkREREQRgsGOiIiIKEJwVywRERFRhOAWOyIiIqIIwWBHREREFCEY7IiIiIgiRNgPUCyKIoqKihAdHd3swtxEREREkUCSJNTV1SE1NRUqVfBtcmEf7IqKipCenq50G0RERETtqqCgAGlpaUFrwj7YRUdHA5B/2MMv10NEREQU7mpra5Genh7IPMGEfbBr3P0aExPDYEdEREQR61gOOePJE0REREQRgsGOiIiIKEIw2BERERFFCAY7IiIiogjBYEdEREQUIRjsiIiIiCIEgx0RERFRhGCwIyIiIooQDHZEREREEYLBjoiIiChCMNgRERERRQgGOyIiIqIIwWBHREREFCEY7IiIiIgihEbpBsKJw+Fo9TG1Wg2DwXBMtSqVCkaj8Yjaeo8fpbUuVDrcEEUJEgBBUEFvMEK+B7icTkiSBI1ahTiTFvFmPSxGLVQqAYIgwGQyBdbrbKhtyeG19fX1EEWx1Z7NZvNx1bpcLvj9/pDUmkwmCIIAAHC73fD5fCGpNRqNUKnkv3E8Hg+8Xm9Iag0GA9RqdZtrvV4vPB5Pq7V6vR4ajabNtT6fD263u9VanU4HrVbb5lq/3w+Xy9VqrVarhU6na3OtKIqor68PSa1Go4FerwcASJIEp9MZktq2/LsPxWfEsdS25d89PyP4GcHPCNmJfEY0/h51GlKYs9lsEgDJZrO1+2sBaHUaPXp0s1qTydRq7cmnny09+9ke6V9v/SqNe3mjpDVbWq3VJfeSMu9bH5jUMd1arTUnZUnXL/tJmvL6Funh936XkjNParU2MzOzWb9DhgxptdZqtTarHTFiRKu1JpOpWe3o0aODvm9NjRkzJmit3W4P1I4fPz5obVlZWaB2ypQpQWtzcnICtTNnzgxau2PHjkDt7Nmzg9Zu2rQpUDt//vygtRs2bAjULly4MGjt+vXrA7UrVqwIWrtu3bpA7bp164LWrlixIlC7fv36oLULFy4M1G7YsCFo7fz58wO1mzZtClo7e/bsQO2OHTuC1s6cOTNQm5OTE7R2ypQpgdqysrKgtePHjw/U2u32oLVjxoxp9jscrLYtnxEjRoxoVmu1WlutHTJkSLPazMzMVmv79evXrLZfv36t1vIz4tDEzwh56gyfEQ8+9KBk99ilWnet9PO2n4PW3nn3nVKxvVgqqiuSftrxU9DamybeJO2r2iftrdor/bQ3eO0Vf79C2l62XdpWuk368cCP0tbSrZLH55HaU1uyDrfYKeBAuR3/+XJf4L5flFqt1WtVOKlbFBr/HihVC2jt71WPT8RPByoD9yvtrf+VVufy4dOdJejVLQqZCeZW64iIIpEoifD6vfBJPnj9rW8pA4ADtgNwVbngF/0oc5YFrf2p6CcUmAsgiiL2VO0JWrv+z/X43fg7REnEz8U/B61944838JPxJ/glP37OCV770m8v4UvzlxAlETt37Qxau2j7IqyPWg9REpG7NTdo7X+2/gdvW94GAJT8UhK09vktz2PtmrWQIKF6Z/VRa1e/uhqSJMFxoPWt0wDw4u8v4n9v/A8A4DrY+pY9AHhrz1v47r/fAQA85a1/HwLAxzkf49f//QoA8NW2vvUWAL4t/BY3f3QzAEB0i1DpVdhw3QZYjdagz+sogiS1sh2+Ay1evBhPPfUUiouL0b9/fzz//PM499xzj+m5tbW1sFgssNlsiImJadc+j7Y7ZG+FGz8fqMSWvGps3l+ESvuRHxaZCUZkWaOQlhiHFIsByTEGWLR+JMcY0M1iQLRe02yz7tF2s3h8IqqdHlQ5PKhxeuEQ1ai0y/dzS6uQW+FAXpUTNudhvQiASivvFjJq1TgpXofeSWb0TY5Gn+Ro9E6KRrRBGyjnbpa213I3iyzcd7NwV6ystX/3kiTBJ/mgM+jg9rvh8XtQ66iF2+uGR/TA4/fAJ/rgET3wil54RS80eg28ohcevweOegc83kOPNQYtn+iDV/RC0AuB2656F7y+Q4/5JB98og9+0Q+v6AV0gF/yy497vPB4PfCLfvgkucYn+uCDT/45dAh8RoheEWj9Iw2CVoCgaqj1iWj1r+sOqpV8EiR/61/dgkaAoD6OWr8EyRekVi1A0BxHrShB8gaJGmpApVGFrFYlyI+r1CqotWoIEORtbd5DjwsN/0GQ7wdqBblW8khH1AgNm1jUGjVU2kP3VYIKr41+DfGG+Nb7PkFtyTqKB7u1a9di7NixWLx4Mc455xwsW7YML730Enbt2oWMjIyjPr8jg11L9pTU4f3tB/G/X4tQWN38i0WnUeGU7hYMzozD4Mw4nJ4ZB2uUvsN7BIAapwc5FQ7kVDiQW+FATqUTORV27Cu1w+1r+RMtPd6Ik5NjMKC7BQO7WzCguwWJ0cr0T0TBSZIEn+hDvb8eLp8Lbp8bLr8Lbr8bLp9Lvn3YMrffHZhcPhc8fg9c/uZzt98Nr98bqPP4PfCInkO3/R40HgMcCdSCGhqVBipBBY2ggVqlDtxWqVRQC2p5UslzlXBoWePjgXpB1WxZ47xxam154zJBEJrPIcj94FBd42ONwePw9QiCHE4a19NYF3hOK4833lahoaZJD02fd/jrNO2jpdqm627pdmBZk9cN3G5Sc/hzGl8nUoVVsDvrrLNw+umnY8mSJYFlJ598Mq666irMmzfvqM9XYotdYZUTH/5ejA9/K8a+MjsAQFCpEGUy4pyTrBiSFYd+iXr0S42BXqM+Yj2d6a9xvyghp8KB7Tml2F1sw+6SOuwpqUNpbfOtNSqdvKUhKUaPkxP1ODk5Gv1TY9A/1YJuMYZmtdxiJ+MWOxm32MlUKhXUOjXqffWo99WjorYCLp8rEMScPidcfpccwkQXfIJPftxXjzpXXSCc1fvq5SDmdwWmxmDml4Js8ukgGkEDrVoLrUoLnUonz9U6aFQaaFXycoPWEHhcBTkENX1co9LI61FpYTKYAo/BD7ledai+MYhpBA2ijFHQqOXHRJ+8FS7wuEoj324Ia9GmaGjU8jK/1w9JlAIB4nCt/buXJAmSBIiSHG1FSYLBYIQgqCBKEtweDzweT8PBgpC3BEGCKMnP1RsMUKnUkCDB4/bA4/UGaiQJjQcZAgB0+qafJ154PO6GHuQem9ZqG/7dS5L8GeFtobZxrtPpoW74nmr67/5QrRS4r9U2/4xo7KFR068bjVbbpNYPj9vV7LUPPUeCRqOFtslnhNvtOqymyXo1Wmh18npFUYTbVX9ETePrNF2vKIpwHfYZ0fT7Ua3WQNvk372r3tniegFApVZDd9jJE0My46HTtN9AI23JOooeY+fxeLBlyxbcf//9zZZfdNFF+PHHHxXqqmU/rluABc8/hLf31UOsFyG6xGabt4ecez6+/eIzGHXyPxCz2dzqF8KIESPw9ddfB+5nZWWhoqKixdohQ4Zg8+bNgfv9+vVDXl5ei7X9+vXDzp2Hjqc444wzsGvXrhZrMzMzkZubC7VKwEndonDDpaPwyy+/tFhrjI7F0Nnv4kCFA6W1bvy69B9wF+xosVZvNOJAUSVSLAYIgoBrrrkGH330UYu1QPN/WGPHjsV///vfVmvtdnsgCN5xxx1YtWpVq7VlZWVITEwEANxzzz1YvHhxq7U5OTnIysoCAMyaNQtPP/10q7U7duxA//79AQCPP/44HnnkkVZrN23ahDPOOAMA8J///Af/+te/Wq3dsGEDRo4cCQBYvnw5pk2b1mrt+vXrcemllwIAXn/9ddxyyy2t1q5btw7XXnstAODdd9/Fdddd12rtihUrMGHCBADAp59+issuu6zV2oULF2Lq1KkAgO+++w6jRo1qtXb+/Pn45z//CQDYunUrzjzzzFZrZ8+ejTlz5gAA/vjjDwwYMKDV2pkzZ+Kpp54CAOTn56NHjx6t1k6ZMgWLFi0CAFRUVKBbt26t1o4fPx4rVqyAy+9Cha0CJ518ElQG1aFJf2h+6pBTcdV1V8Hpc8LpdWLl6ysh6AWoDepmdSq9CmqjusMGmJJECaJHhORpmHvluUFjwNlDzoZBbYBercenH36K6vJqiF65prFO8kqwRMXi4dlPQC1ooRa0mDfrQRzYs0/eteeV5Oc03I4yW/Dqx79AgAZ+UcCs28Zg55aWP8N1BiOeXv8rfKIEv0/CSw/fjj2bv231Z5mxZhv8ogS/KOGT//wLBzZ/2WrtFc99AUFrgCS5sfW1uSjc9EmrtWfMehsaswWiBBz4339QvvGDVmuzp62A2pIESZJQ/uXLsG18p9XalFsXQZeYCQCo+f512H54s9Xa5HHPQp/SGwBg2/g2ar5e0Wpt0g2Pw5BxCgCgbut6VH2+tNXaxDGzYeopf/bYf/8ClR8932qt9cr7Ye47HADg2P09Kt5/otXahNEzEDXwLwAA55+bUf7f1j//4i+cjOjT5c8QV/5vKH3zwVZrY0feAstZ1wAA3MV7UfLqPa3WWs65AbHDbwIAeMrzUPzK1FZrY868GnGjbgUA+GylOLh0Yqu1UYMuRcJFdwIA/E4bChfc1GqtecAFsF76DwCA6HFBpTNg86y/dJo9WooGu4qKCvj9fiQlJTVbnpSUhJKSlg/MdLvdzbYk1NbWtmuPAHBw7zbo5y7G/7lTMT0D2JcqYG93AXuSJOyzSHD6/PDrqnHnVxNhNVqRYk5BzMgYqIvV8FZ64a30wu9Q/i/p42XWa/DlvSNhd/uwq6gWYz81Y29By7Uen4hhT3yFBLMO/btbcKDc3rHNUtelhhyoDCqojHKoUhvl+yXWEqzetRoOrwMVtRVIGZsSeKzppDaosT16Owa9Niiw9av3U71bfclylOPF318M3LecZTm2VgUtPHYXRJcI0SP/oSi6GyaPiChTPPqf8ReooIMKenzz5stw1zoCj0tuOYCJHhEWaxbOu/1JSKIWfr8Wnzx8M5zlLX9+GhIzoZl2HXx+CV6/iJyV/4W3oqjF2poYNeb0btxbIKF4sxOekpb/WK31eHHXG4f+qCwpa/3fvdcv4pnP9wbul1W2vkUUAN7ddjBwu7w2+MHy2/JrAnsXqg8/rvgwxbUuqH3yF3G9J8gBdgDcPhGahkNW2mMflyDIUzBatQqGhuO6XOrgfyGYtCpYjA1btXRqVAapjTFoEBelAyBAbdCg5U0MDbVGDRJj5PesxqhFeZBai0GL5IY/8G1VepQGqzVqkRYn/67V2fUIdlqGxahFRrwJggA4vAYUH6U22ypvDKhXmXAwaK0GPRPNEAQBHrsPhUFqYwwa9E6KAgD43Wqo9UaoVZ1nN7Ciu2KLiorQvXt3/Pjjjxg6dGhg+dy5c/Haa69h9+7dRzxnzpw5LW4lac9dsc6SYnxzw/lIKwU0h/379wtAXtKhsLc9W0CdqYVN+Wojks3JSDImITMmE70TeqOHpQd6xvaE3t96ylf6wOhGrZ084XD7sKekFruKarGruA47i2zIsfkDZ/pKPg+khlqzTo1eSVHolRSNvsnR6NUtCoOyU2AxyR9A3BXb9XbFSpIEu8eNalcdalx1cPjr4fC7UOupg81VixpHNew+R8MWMQfq/U64/E7U+x1wi/XwSPVw++W5Xwp+1tvxkkQdIOog+XWQms5FHSAaABggiXpA1EH0qiD5dYCoDdQ0Ph+SHlBFybehhugJElIEASrtoc+FNtV6XUfu7wrUHjpp6mi1Wo0And4EjUqAWi1A5fdCLUjQqAVoVAI0KhVUgiDfV6tgMJqgVsmPST63/DyVALUgHLrdMBlNZqhVgEalgujzQAUxUKsSmteaTGaoGtbr87ohSBLUKjSrUwnyc40mEzRqFVQC4PN5IPn9DesEVI11qsbPP3PguT6PB6LkgxoCBEFet9D4PEGA0WSERq2GIAA+rxc+rzfweGAOAYIKMBmNUKvVUAmA1+OFz+eVg1vDugU0rl/+POFnRPgfrtERx/aFzTF2Ho8HJpMJb731Fv72t78Flt99993Yvn07vvnmmyOe09IWu/T09HY/xs7123vQvTUeLnssHH3+CfvOP+D+9TdIZc3/vvEbddh9xSn4aXg8DrpLUewoRpWrKui6Y/WxyLZkIzs2G9mWbPS09ETv+N6d5tTptnJ5/dhdUocdB23YWWTDjoO12FNSB4+/5eCYajEEzsTNTjSjhzUKPaxmWKN0EX0wbGfn84twePxwuH1wenywu/1wun1wevxwev2o9/jgdPtg8zhgc9tQ565DnbcODm8d6v0OuPx2uPxOeCQHvKITXjjhhxMi6iEJ9ZBULkDlgqAKPrRAW0mipiFo6SE1TBD1kPz6lpeL+obwdeRyiFocy/7TxkCjVaugVgnQquXgo1HLyzQqOfxom4SixkCkVQmH3Va1/vzD1qNu9TlNXqNx2VEeU6sP/QyNgYWIOoewCXaAfPLE4MGDmx0D1a9fP1x55ZWd6+QJUQSWDAPK/wDOfwg4Tz5myFtcjPrt21G//Vc4fvwR7n3y+HS6zEx0u/8+RI0cCbffjRJHCYodxThoP4hcWy7+tP2JHFsODtpb3zgcb4hHn7g+6B3XG33i5Xm2JRtatbbV53RWXr+I3AoH/iipw56SWuwursPukjocrGn9L6hogwbZVjN6WBvCXqIZ6XFGdI81whqlh6oTbfruTNw+P+pcvobJizqXD7X1XtS5fbC7fLC75amu8bbLi1p3PWrdNjj8tXD67HD56+CDA4K6vmFyynOVE4LaBUFdD6jq5dtC8N1Yx0ry6wDJAEE0QJCMUMEAtWSEGkZoBHnSCkboBBN0KiN0KiP0ajMMahP0KhOMGhNMGjN0ai20mqaBp0mgOiwA6dQqaNSqZluhtA0hqzHkaJo+t2EeCEMNgY6/i0TUnsIq2DUOd7J06VIMHToUy5cvx4svvoidO3ciMzPzqM/v0OFOfnsLeOc2wBgP/GMHoGs+sK8kirC99z7Knn0W/oaTIcznnoukB+6HPju7xVU6vU7k1ubigO0ADtQcQI4tB/tr9iOvNq/FIQQ0Kg16WnoGgl7f+L7oE9cHsYbYkP+4HcFW78XeUjnk/Vlmx4EKB3Iq7Cisrg96LItWLSApxoBUixEpsQakxhqRajEgxWJEQpQOcSYd4sw6xBg0YbXlweMTUefyBoJXbUMwawxptfVNwlrgMS9qm8w9Pg8EjROC2g5B7YCgcchztbPJ3NlQ07BMFfx4pKMRoIZeiJKDlioKRk0UTBozTJoomLXRiNZFIUYXjRh9DGL10bAY5HmcIQaxxmjEG2Jg0IbX/ysioo4SVsEOkAconj9/PoqLizFgwAA899xzOO+8847puR0a7Pw+YOEQoDoHuGguMKzlMxf9djsqly5F5apXAa8X0GgQf9NNsE6bCnV09DG9VL2vHvur92Nv9V7sqd6DPVV7sLd6L+zelg9KTjIloU98H/SJ6xOYZ8RkBAZqDDcurx/5VU4cKJfH3jtQbkdOhQMHa+pRWutCkIt1NKNWCYgzaRFr0iHepEOsSYs4kw5GnVqetPJkaHLbqFNBr1FDABrGUTp0TIy8YUY+pkaUJHgbDkKXpyNv13v88q5LT8Ouy8B9+Xa9199sy1rLYwqKDUHMLoc1jQOCuq5hboegsUMVCHB2COrWj3sJRgUVzNoYROtiEKOLQazegnhDLOKMsYjRxcCityCm4bEYfQyitXJQi9ZFw6A2MJQREbWTsAt2J6LDByjesgr4YDoQlQzc/SugNbRa6snNRekTT8LeMLSJOj4eif+Ygdirr4agPnJ8u6ORJAlFjiLsrdqL3dW7sbdKDn0FdS2fomrUGJFtycZJsSehV1wv9IztiZNiT0KSKSmsv4R9fhFldW4U1dSjyOZCcU09im0uFDXMqxweVDs9cHo685nIEqCqh0pTB0FTJ4c2TR0EtR0qTR00OjtUDctEwQ4IbftnqhbUsOgtiDfEI84Qhzh9HGL1sYg1xMpzfSwsektgucVgQbQ2Oqx/L4iIIhWDXXvyeYAXTgNqDwKXPguc0fq4OI3s332H0sfnwZOTAwAwDxuG7s8/B3WI+rV77NhXsw+7q3ZjT5W8dW9fzT64/S1vuYnWRuOkuJNwUuxJ6GHpgcyYTGTGZCI1KlUeBDRCuLx+1Di9qHZ6UO3woMrpQbXTC5vTg3qvH/UeEfVeP1xev7xlzeuHq2ErmtvnDwwQKjbcaDoIqSTJwxPoGo7h0mlUDQez+6HS2CGpayEKtYCmDqKqFn6hFl7UwC3ZUC9Ww+GrgYi2nTQQq49FgiEB8cZ4xBvi5duG+MD9eEM8YvWxiDfEI1oXHbZba4mIqDkGu/b281Lgk/uA2Azgrq3AMZzMIHm9qHr9dZT/5wVI9fXQZWcjfcli6I7hOMLj4RN9KKwrxP6a/dhXsw/7q/cHjt1rbXR6taBG96jugaCXEZOBzOhMpESlINmcDKPG2OLzIp3H70GVqwqVrkpU1cvzyvpKVLoqUVFfgcp6eV5RX4FaT9vGVYzRxSDRmIgEYwISjAmwGq1IMDTMm9yPM8RBo1J02EkiIlIIg1178ziB5wcCzgrgqqXAaTcc81Ndu3ahYMpU+EpKoLZY0H3BCzAHGYk/1Dx+T+AEjcagl1ebh/zafLj8wQf/tOgtSDIlIdmcjGRTMpLM8u0kU1Jg916sIRZ6decYfbslPtEHh9eBGncNql3VsLltqHHXBKbGZVWuqkCYq/PUtek1tCrtoYBmsiLRmAir0dpsagxzOrWunX5SIiKKFAx2HeG7Z4EvHwGsvYEpGwHVse/28paVoXDqNLh+/x3QapEyZzZir7mmHZs9OlESUeYsQ35tPvLq8pBny0NeXR4KagtQ7CiG0xd8dPimjBojLHoLLDpL4FiuGH0MDGoDDBr5UkYGtQF6jb7ZMo1KAwGtH+MlQb7IeeOFyRsvQt507va74fA6YPfaYffY5XmT2/W+1odXCUYjaOTdn8aEZvOmQa1xK1uMLobHqhERUcgw2HUEVy3w/ADAZQOuXQX0v6pNTxddLhQ98ADqPpavZRh/663odu89x3VSRXuTJAl2rx0ljhKUOEpQ6ixtdrvMWYYadw1sblunuAj5sTBrzYe2Mh52UkHj/QRDgjwZE3jMGhERKYbBrqN8NRf4dj6QfApwx7dHv9jfYSRJQsXCRahouDh51KhRSH3qKaijzEd5ZufUGAAbQ17j7k2b24ZaTy3cPnmLWr2vPrB1zeVzBeY+KfjJBJIkQavWQq/WQ6fWQa/SH7rdZB6li4JZa0aUNkqedIfNtVFhOcgzERF1TQx2HcVZBTw3APA6gBvfAnpfdFyrsX34IYofeBCSxwN9nz5IX7wI2u7dQ9wsERERhaO2ZB3uWzoRpnjgjFvl2989jaCXSgjCcumlyHztVaitVrj37EHOdX9H/c6dIWyUiIiIugIGuxM1dBqg1gMFG4Hc7497NcZTT0WPdWuh79sX/spK5N86Ea5du0LYKBEREUU6BrsTFZ0MnD5Wvv3d0ye0Km1qKjJXvwbjaadBtNmQf8utcO3eHYImiYiIqCvgMXZt4HA4Wn6gpgDqpcNgUPuB274C0ga3XgtApVLBaDw02O/htWKdHWVTp8KzcydUsbHIXLUShj59AABOpxOt/S8TBAEmkylwvy219fX1EMWWrlMqM5vNx1Xrcrng97d+pmxbak0mU2AYEbfbDZ+v9ZMt2lJrNBqhahiuxuPxwOv1hqTWYDBA3XCWc1tqvV4vPB5Pq7V6vR4ajabNtT6fD25369eR1el00Gq1ba71+/1wuVofA1Gr1UKn07W5VhRF1Ne3PjxNW2o1Gg30enl8RUmS4HS2PnxPW2rVajUMhkOXFQz2774ttUf7jAhWy88IfkbwM6LttSfyGdERw1u1KetIYc5ms0kAJJvN1u6vBfkKUy1OowdnSNLsGElafa0kSZJkMplarR0xYkSz9Vqt1iNqolUqaU1GprSrT19pz9lDpfo9eyRJkqTMzMxW19uvX79m6+3Xr1+rtZmZmc1qhwwZ0mqt1WptVjtixIhWa00mU7Pa0aNHB33fmhozZkzQWrvdHqgdP3580NqysrJA7ZQpU4LW5uTkBGpnzpwZtHbHjh2B2tmzZwet3bRpU6B2/vz5QWs3bNgQqF24cGHQ2vXr1wdqV6xYEbR23bp1gdp169YFrV2xYkWgdv369UFrFy5cGKjdsGFD0Nr58+cHajdt2hS0dvbs2YHaHTt2BK2dOXNmoDYnJydo7ZQpUwK1ZWVlQWvHjx8fqLXb7UFrx4wZ0+x3OFjt6NGjm9We6GdE4zRkyJBmtfyMkPEzQsbPCFl7fkZ0hLZkHe6KDZW4HoCgBvZ9Cvyx/oRXVyeKmFRYgByVCv7qauRPuAXufftC0CgRERFFKu6KbYOj7mb5/gng++eA6BQ4JnwFGCwt1rZlN4tkt6P8zilw7doFdUICui1bCm2PHi3WcjfL8dVyN4uMu1naXstdsYfwM6LttfyMkIX7Z0Rn2xXLYBdK3npgyTCg6gAw5FbgsudCslp/TQ3ybrkV7j/+gNpqRearq6DPzg7JuomIiKgJUQR8rsMmd8Pcc9j9hvkp1wGa9rtOOoOdknK+BVZdLt++5RMgc2hIVutr3B27Zw/UiVZkrnoV+uyWt9wRERFFHJ9HviCAxwl4G6f6Q3NPC8t89Q33XU1u18thzOuUg1ngfkOI87e+JbJVM/cDUYmh/5kbMNgp7f1pwLbXAGtv4I7vAK3h6M85Br7qauSPnwD33r3QJCYi843XoUtPD8m6iYiIQkKS5MDksQPuukNzt12+7XE0TC3ddzSfvM6GwOYAxOCXnWwXghrQGuWtcRpD6/MrFwPmhHZrg8FOafXVwMIzAUcZcN6/gPNnhWzVvqoq5I8fD/e+/dBmZiDrjTegSWi/XyYiIupC/D7AXQu4bPIUuF3bEM5qG6a6JsvqDi1rDHLtGcJUGkBrBnQmOXRpm84PX2YANMaG20Y5hGlNTZY3zDX6Jo83BjkjoNa038/RBgx2ncHO94C3xsu/gHd8ByT1C9mqvaVlyLvhBniLimAYMACZq1ZC1eQAYyIi6sJ8bnkDQ7OppiGg1bRyu2Hytn6iznHRmgF9FKCLaphHN8zNDVMLt7WmhtsNQS2wzCSvT6MLbY9hgMGuM5AkYM1NwJ4Pge5DgImfASp1yFbvPpCDvBtvhL+mBuZzzkH6ksUQdF3vl52IKGJJkhy2nJVyOHNWNkxVh27XV8nhrGmI87Z+Rucx05oBQ4w8uoM+Rr6tjwH00c1v66MPux3dJMRFhfR7rytjsOssaovkXbKeOuCS+cBZd4R09fW//oq8CbdAqq9HzOWXI/XJJyCoODQhEVGnJEnyVjJ7OeBomJwVgKOyYV7RsKxSvu2sBKTWh3cJSlABxjh5MsQCxlg5pAVuN9xvutwQI8/10YBaG5IfmUKDwa4z2fwS8OG98l8/UzcCsaE92cH+7bcomDIV8PkQf8stSLrvXyFdPxERBSFJ8vFldaWA/bDJUQHYy+TjrRvDnNj6WHWt0kUBpnjAGA+YEuTbpoSG+/GHAlzTSR8D8A/9iMFg15mIIrDiEqDgZ6DXRcCN64AQD2ZY8957KL7/AQBAt3/9Cwm33hLS9RMRdTmNga22GKhrMtUWA/aS5kHO1/qAui3SWwCzFTAnNsytgKnpPKHJ/YR2HR+NwkNbsk7nON0jkqlUwBUvAEuHA/s+A3a8DQwcE9KXiL3qKvgrKlD29DMomz8fGmsCLFdcEdLXICKKGKIo7+asPdgwFQG2QnleW3QoxLXlWDV9DBCVBEQnA1Hd5NtmK2DuJt83J8pzkzVkQ2ARtYTBriMk9gHOnQl8/Tjw8X1Az/PlzechFD9xInzl5aha9SqKHpwFdVw8os4dHtLXICIKCx6HHNRsBUBNwaHbtsJDQc7f+mW2mjFYgOhUObBFpwAxKUBUMhCdJM8bQ5zOdPR1EXUA7ortKD4PsOw8oPwPoM+lwN9fC/nZQpIoouif/0Lthx9CMJmQuWoljAMHhvQ1iIgU53EANflAdR5Qk9dwO/dQeHNWHsNKBDmUxXQHYlIBS5o8j+kuB7jGIMfARp0Aj7HrrAq3ACsulv9SHDwBuOz5kB9vJ3k8KJg8GY4ff4I6Ph5Za96ELiMjpK9BRNSuRFHeFVqdI197uypHDm41+XKQc5QffR26aPlkNUu6HNoab8d0Byzd5a1tXXA8NApPDHad2a73gXXjAUjAiPuAUQ+G/CX8dgfyx42Da9cu6DIzkbnmTWji4kL+OkREx030y1vYKvcDlQcaQlxDkKvJO/oJCXoLEJcBxGYCcVnyvGmQM8Z2xE9B1CEY7Dq7zS8DH94j3x79NHDmpJC/hLesDLnXXw9fUTGMp5+OjBWvQKXnmVVE1IEkSR7yo3I/ULmvYf6nPK86EPw4N0ENxGYA8T2A+OxD4S0uU15u5B+r1HUw2IWDDfOAb54AIADXrgD6/y3kL+Hetw+5N94Esa4O0ZdcjO7PPMMBjIko9ES/vJu0Yi9Qvgeo2AOU75XnLlvrz1Pr5NAW37MhwPUA4hqCnCW901ynk0hpHO4kHIy8Xx608pdXgHdulweazB4R0pfQ9+qFtAUvIH/S7aj7+BOUd++ObjNnhvQ1iKgLEUV5N2nZrobpDznAVe4LsutUkHeRJpzUZOopzy3pvOQUUYhxi52SRD/w1njgjw/kA31v+RBIOTXkL9N0AOPkObMRd/31IX8NIoogkiQPvNsY3kobglz57tbHdlPr5bCW2Buw9jk0TziJ47YRnSDuig0nXhfw+hgg9zt5AMuJn8m7IUKsfNEiVCxYCKhUSFu8CNEjR4b8NYgoDPl98jFvJb8Dpb/L85LfWz/zVK2XQ1u3/kC3vkBiX8DaWz4GjlvfiNoFg124cdmAFZfKH6pxWcDEz+XxlUJIkiQUPzgLtnfflce4e/VVGAf0D+lrEFEn53UBpTuBoq2HAlzZrpZ3owoq+Y/Mbv0appOBpP7yMXA89o2oQzHYhaO6EuDli+TjV5IHAuM/CPlZX5LXi4I77pDHuEu0oseaNdB27x7S1yCiTsLnbghx24Di7fK87A9A9B1ZqzUDyQOApAHy50/yKXKQ4+C8RJ0Cg124qvxTDnfOCvmg4r8tBbJCe1kwf10d8m66Ge69e6HvdRIyX38d6nB/34i6OlGUd6ce/AUo3Awc3CIfFyd6j6w1JQCpg+TwlnKKPI/rIV/Xmog6JQa7cFbyO7B2rDxYJwRg+Axg5IMhHSHdW1yM3L9fD19ZGUxnn42M5csg6DgCO1HYcFYBB7fKIa5wsxzoWhpWxBgnh7jUQUDKafLckhbyK94QUftisAt37jrg4/uB7avl+ymnAle/JB+wHCKuP/5A3k03Q3Q6YbnyCqQ88QQEftgTdT6SJG/Nz/8JyP8ZKPhZ3jp3OI0RSD0NSBsCdB8MpJ4uD+TLf9dEYY/BLlLseh/433TAVSN/aP91LjDk1pB9UNu/+w4Fk+8E/H5Yp0xB4vS7QrLeDuH3yeMA1pXIwzI0zu1lgLdeHtHe7wZ8LcwlP6A1ATozoIsC9FENtxvu66IAg6Xh+pIZ8rUlOVwDdRS/Dyj5VQ5xjWGupTNU43sCaWfIQS7tDPnEBrW24/slonbHYBdJaouA9+4EDnwt3+99CXDFAiAqMSSrr163DiX/NxsAkDJ3LmKvuTok6z1hgUsR7ZNHs6/YJ0+1RYC9RH4MHfira+522AXFM+Qv1uQBQFQSt4rQ8fN55GPicr8Dcr8HCn8BvI7mNWq9vBUu42x5SjsDMMUr0y8RdTgGu0gjisDGJcAXc+QtUeZE4IqFQO+/hiRQlD37HCqXLwc0GmQsXwbzsGEn3vOxkiQ5rJX8DpT/cSjAVeyVt1QGI6jlYWGikoDoZHkelSSfyafWy8clqvWARi9fuqhxrlLLW/XcdsDTODnkubvhtrMSsBXKFylvbUDWRiZr8zMKkwbI43qF8LhIiiB+n3yGau63QM53QMHGI3/HDLGHQlzGUPnYOA2v9UzUVTHYRaqSHcA7k+RxpwAg8WTg9LHAKdcD5oTjXq0kiij6132oXb8eqqgoZL7+Ogx9Qnc8X4DfJwe2kt+Bkt8OjaNVX9XKExouRWTtDST0AqwnyRcBbwxypoT2HxBVkuQD1W35ctCrKZDDXuN1MSv3A5J45PNUWnng1u6ny1/MmUPl3rllr+sRRaB0h7zVPedbefeqx968xpQgnwGfdS6QeY78u8OzVImoAYNdJPO6gK8eAza/dGhQUZUW6DsaGDQO6DnquMKO6PGgYOJtcG7eDE1yMrLWroE2Ken4+3TVymNoNYa40h3y8At+95G1glr+IkvqJ4c4ay85yCX0BLTG4++hI3ic8pbGkh3yz9g4d9ceWRud0rAFZpg8T+rPkfojVW0xcGAD8OdXcqA7/Bg5Y5wc4HqcJ4c5BjkiCoLBriuorwF2/BfY+po8+GijmDTgtBuBQTfJV7FoA7/NhtwbboTnwAHo+/ZF5urXoI6KCv4kUQRqC+VA0/SSRNW5LdfrouXdlsmnNAyEOlD+UoukkxMkSd6iV/IbULBJ3kJTtP3IMcX0MUD6mfKXe8/z5Us08cs9PHmcQN4PwJ8NYa78j+aPa83yFrnsEXKQSxrA/9dEdMwY7Lqakt/lgPfb2ubHpaWcJm/9is9uPpkSWt0l6CksRO7fr4e/shLm4cORvmQxBEGSdz9W5QBVB+Qx9qpy5Hl1bsuXIwLkkJnc5LizlFOA2Kyu+YXmccqXccr7SQ56BZsAT13zGnM3eYtr9ih5Hp2sTK90bCr/BPZ9Jk+5Pxy2NVqQj4vrOUoO7Wln8phLIjpuDHZdldcF7F4PbHvt0Fm0LdFbgPge8qQ1ycHM5w7M6wtqkbe2HJIPiO0LJJ9WAgEtHEfWqPF4ssYQ1xjkeNZe60S/vKu6cStP7ndHHkDfrf+hYJB5TmRt1QxH3no5wO37DNj/ufxHTlOW9EP/v3qM4O8/EYUMgx3JB/kXbZO/fAJTjrzb9BjUHdSj8Pt4QBKQOLAW1tNE+bJD8T3kXbzxPeStf3E95C80XhT8xPjc8la8P7+Sj80q2o5mw7loTfKWvN4XAb3+CsSkKNVp12IrBPZ+Auz9TD7xwVd/6DGVVj4pptdF8mTtzZNjiKhdMNhR67z18u7TxqAn+gCNQR5K4bB51Uc/onTxGwCAlLn/Ruw11yjbe1fiqARyvpaD3v4vgbri5o8nnwL0vlieUgd1zd3b7UEU5T+I9n4M7PlEPma0qehUoNeFcpDLHgHoo5Xpk4i6FAY7CpmyZ55B5YsvAWo10hcvQtSIEUq31PVIknwc5d5PgX2fygPYNt2aZ06Ug0afS+TdgDqzYq2GJY9DPnRhz8fyblZ7aZMHBfkEl95/lbeUJvXnVjki6nAMdhQykiSh+P77YXv/fxCMRmSuWgnjKaco3VbXZi8H9n8h7yL886vmQ6uo9UD2SHn4m94X8wSM1tSVyO/f7o+AnG+anwCkiwZOOl++ykuvi05ojEgiolBgsKOQkrxeFEy+E44ffoA6Lg5Za96ELjNT6bYIAPxe+SzbPZ8Aez48cpiZ7kPkLXl9RgPdTu66W5skCSj7Q36P9nwsX8KrqdgMOcj1uRjIHM4zWImoU2Gwo5Dz2x3IHz8erp07oU1PR9abb0BjtSrdFjUVCC8fydPh4SUuC+h7mRzyMs6O/MGR/V4g70c5yO35CKjJa/54IPReAnTr13VDLxF1egx21C58FRXIvf4GeAsLYejfH5mvroLKzOO5Oq3aYnl3456PgAPfNB9nzZQg76rtM7rhuDyTcn2GkqMC2Pd5y7upNQZ5N3WfS7ibmojCCoMdtRtPbi5yb7gR/urqQwMYa7VKt0VH47YDf34pH1O295PmA1lrjPL4a30vlQOPOYy2xEqSfAm3xiFJCjej2YklJmtDgL1E/hl5YgkRhSEGO2pX9b/9hrzxEyDV18Ny5ZVIeWIeBO7GCh9+H5D/oxzy9nwoX/4sQADSz5JPvugzWr5ySWfjqJQHdD7wtXwWa+3B5o8nD2wyFMzpHAqGiMIegx21O/s336BgylTA70fCpEnodu89SrdEx6Nxi1djyCv+tfnjCb0Ohby0M5Q5Ls9tl08QOfC1fAZryQ402yqnMcq7WHv/VT6L1dK943skImpHDHbUIWrefgfFs2YBAJIefADx48Yp3BGdMFvhoZMNcr4DRO+hxxp3a2aPADKHAZa09unBUSkPDJz3o3xs4MFf5IG0m0o8We6j5wVAj3MBrbF9eiEi6gQY7KjDVCxdivLn/wMASHliHmKvukrZhih0XDb5qhd7PpJ3ebpszR+3ZMiX1MocBmQMk3fbtmWXvN8HVO6XtxiW/C7PS3ceeZUNQB6OpMcIectcj/OAqG4n9KMREYUTBjvqMJIkoeyJJ1C16lVArUbaC/9B9AUXKN0WhVrj0CH7PpPnxb8Ckr95jckqD6OSOki+73PLA//63PI1VpverysGynY3P1O3qbgeQPfTG8LcCHmoFiKiLorBjjqUJIoonvUQbO++C0GnQ/ry5TCffZbSbVF7ctuBwk1A3k9y0Dv4S/OrNxwrrVm+TFdSfyB5AJA0EEjqx2uwEhE1wWBHHU7y+VA4YwbsX3wJlcmEjFUrYRw4UOm2qKP43EDRdvls27LdgFojn9Sg0cvjx2kMh25rDYAhVg5zcT141ioR0VGETbDLyspCXl7z0eDvu+8+PPHEE8e8Dga7zkN0u1Fwx2Q4f/4Z6thYZK5+DfqTTlK6LSIiorDWlqyj+J/Kjz76KIqLiwPTQw89pHRLdJxUej3SFi6EYeBA+GtqkD/xNngKDx79iURERBQSige76OhoJCcnB6aoqCilW6IToI4yI335MuhO6glfaSnyJ94KX0WF0m0RERF1CYoHuyeffBIJCQk47bTTMHfuXHg8nqD1brcbtbW1zSbqXDRxcch4+WVoU1PhzctH/m2T4Of/JyIionanaLC7++67sWbNGmzYsAHTpk3D888/jylTpgR9zrx582CxWAJTenp6B3VLbaFNSkLGilegtlrh3r0bBXdMhuh0Kt0WERFRRAv5yRNz5szBI488ErRm8+bNGDJkyBHL3377bYwZMwYVFRVISEho8blutxtu96Gxr2pra5Gens6TJzop1+7dyBs3HmJtLUxnn430pUugMhiUbouIiChsKHpWbEVFBSqOckxVVlYWDC18uR88eBBpaWn4+eefcdZZxzYOGs+K7fzqt29H/q0TITqdMA8fjrRFC6HS65Vui4iIKCy0JetoQv3iVqsVVqv1uJ67bds2AEBKSkooWyKFGU87DenLlyF/0u1wfP89Dt49A2kv/AeCTqd0a0RERBFFsWPsfvrpJzz33HPYvn07cnJysG7dOtxxxx244oorkJGRoVRb1E5MQ4YgfcliCHo97F9/jYP33gvJ6z36E4mIiOiYKRbs9Ho91q5di5EjR6Jfv374v//7P0yaNAlvvvmmUi1ROzOffTbSFi2CoNWi7vMvUHTffZB8PqXbIiIiihi8pBh1uLqvv0bhXdMBrxeWK69AyuOPQ1CrlW6LiIioUwqrK09Q1xM9ciTSnnsW0Ghge/9/KJ49G5IoKt0WERFR2GOwI0VE/+Uv6P70U4BKBdt/30bJY48hzDceExERKY7BjhQTc/HFSH3yCUAQUPPmGpQ+Po/hjoiI6AQw2JGiLJdfjpS5cwEA1a+9htK5jzPcERERHScGO1Jc7NV/Q8q/HwMEAdWrV6PkkUd4zB0REdFxYLCjTiF2zBikPP64vFt2zVqU8IQKIiKiNmOwo04j9m9XIXX+k4BKhZq3/oviWQ9B8vuVbouIiChsMNhRp2K5/HL5bFm1GrZ330XRAw9wEGMiIqJjxGBHnU7M6NHo/swzgEaD2v99gKJ/8QoVREREx4LBjjqlmIv/irTnnwO0WtR+9BEO3sNryxIRER0Ngx11WtF/+QvS/vMf+dqyn32Gwn/8A5LHo3RbREREnRaDHXVq0eePQtqihRB0Oti/+BKFd02H6HIp3RYREVGnxGBHnV7UeechbfFiCHo97N98g4I7JsNvdyjdFhERUafDYEdhIWr4OUh/cTlUZjOcGzcif+Kt8NtsSrdFRETUqTDYUdgwn3kmMlaugMpigevX35A3bjx8FRVKt0VERNRpMNhRWDEOHIjMV1+F2mqFe88e5N08Ft7iYqXbIiIi6hQY7CjsGPr0Rtbq16BJTYEnNxe5N90ET16e0m0REREpjsGOwpIuKwtZq1dDl5kJX1Excm++Ga69e5Vui4iISFEMdhS2tKmpyHx9NfS9e8NfXoH8seNQ//sOpdsiIiJSDIMdhTWN1YrMV1fBcOop8NtsyJ8wAc7Nm5Vui4iISBEMdhT21LGxyHj5FZjOPBOiw4H82yah7qsNSrdFRETU4RjsKCKoo8xIX74MUaNGQXK7UXjXXbC9/77SbREREXUoBjuKGCqDAWkv/AeWK68A/H4U3Xc/ql59Vem2iIiIOgyDHUUUQatFyrx5iB8/DgBQ+vg8lL/wAiRJUrgzIiKi9sdgRxFHUKnQ7f77kXj3dABAxeIlKH3sMUiiqHBnRERE7YvBjiKSIAiw3nknkmf/HyAIqH7jTRTN/Cckj0fp1oiIiNoNgx1FtLgbbkDq008BGg1qP/oIBVOnQXQ6lW6LiIioXTDYUcSzXHop0pcshmA0wvHdd8ifeBv8NpvSbREREYUcgx11CVHnnouMl1+GKiYG9du2Ie/msfCWlindFhERUUgx2FGXYTp9EDJfew2axES49+1D3o03wpObq3RbREREIcNgR12KoU9vZL75BrSZGfAePIjcm26Ga9cupdsiIiIKCQY76nJ0aWnIev116PudDH9lJfLGjoNj4yal2yIiIjphDHbUJWmsVmSuWhW4vmzBpEmo++ILpdsiIiI6IQx21GWpo6OR/uJyRF/4F0geDwqn342a//5X6baIiIiOG4MddWkqvR7dn3sOljHXAKKI4oceRsWLL/ISZEREFJYY7KjLEzQapDz2GBImTQIAlD/zLMqenM9LkBERUdhhsCOCfAmybvfeg2733QcAqFq5EsUPzoLk8yncGRER0bFjsCNqIuGWCUh5Yh6gVsP23nsovHsGRLdb6baIiIiOCYMd0WFir7oKaQtegKDTwf7llyiYdDv8drvSbRERER0Vgx1RC6LPPx/pL70IldkM56ZNyB8/Ab6qKqXbIiIiCorBjqgV5jPPRMarq6COj4dr507k3XgTvEVFSrdFRETUKgY7oiCM/fsj8/XV0KSmwJObi9wbb4L7zz+VbouIiKhFDHZER6Hv0QNZb7wBXc+e8JWUIO+mm1H/+w6l2yIiIjoCgx3RMdAmJyNz9WswDBwIf00N8sePh+Pnn5Vui4iIqBkGO6JjpImLQ8aKFTANPRui04mCSbej7quvlG6LiIgogMGOqA3UUWakL1uG6AsvhOT1ovCu6bB9sF7ptoiIiAAw2BG1mUqnQ/fnnoXlyisBvx9F//oXqteuU7otIiIiBjui4yFoNEiZ9zjibrwBkCSUzJ6NypdfUbotIiLq4hjsiI6ToFIh6eGHkTBpEgCg7KmnUP7CC5AkSeHOiIioq2KwIzoBgiCg2733IPEf/wAAVCxegtJ58xjuiIhIEQx2RCFgveN2JD30EACg+tXXUPzww5D8foW7IiKirobBjihE4m++CSnz5gEqFWz/fRsHZ86E5PEo3RYREXUhDHZEIRT7t6vQ/bnnAK0WdR9/gsK7pkN0u5Vui4iIuggGO6IQi/nrRUhfvAiCXg/7N9+gcMpUiC6X0m0REVEXwGBH1A6izj0X6cuWQTAa4fjhBxRMvhOi06l0W0REFOEY7Ijaifnss5Dx4nKoTCY4f/4ZBbffAb/doXRbREQUwRjsiNqRacgQpL/8ElRRUXD+8gsKJk2C325Xui0iIopQDHZE7cw0aBAyVrwCVUwM6rdtQ/6tE+GvrVW6LSIiikAMdkQdwDhwIDJWvAK1xQLXb78hf8It8NfUKN0WERFFGAY7og5i7N8fGa+ugjouDq5du5A34Rb4qquVbouIiCJIuwa7uXPnYtiwYTCZTIiNjW2xJj8/H5dffjnMZjOsViumT58ODwd1pQhl6NMHma+ugtpqhXv3buSPGw9fRYXSbRERUYRo12Dn8Xhw7bXX4s4772zxcb/fj0svvRQOhwPff/891qxZg7fffhv33ntve7ZFpCh9r17IfPVVaLp1g3vfPuRNmABfZaXSbRERUQQQpA64WvnKlSsxY8YM1Bx2TNHHH3+Myy67DAUFBUhNTQUArFmzBhMmTEBZWRliYmKOuu7a2lpYLBbYbLZjqifqLDx5ecgbNx6+0lLoe/dGxqqV0MTFKd0WERF1Mm3JOooeY/fTTz9hwIABgVAHAH/961/hdruxZcuWFp/jdrtRW1vbbCIKR7rMTGSsXAFNYiLce/fKZ8vyhAoiIjoBiga7kpISJCUlNVsWFxcHnU6HkpKSFp8zb948WCyWwJSent4RrRK1C32PHshYtVI+5u6PP5A/8TYOhUJERMetzcFuzpw5EAQh6PTLL78c8/oEQThimSRJLS4HgAceeAA2my0wFRQUtPVHIOpU9NnZyFzxCtTx8XDt3In82ybBX1endFtERBSGNG19wrRp03D99dcHrcnKyjqmdSUnJ2Pjxo3NllVXV8Pr9R6xJa+RXq+HXq8/pvUThQt9r17IWPEK8sdPgOu331Aw6Xakv/QS1FFmpVsjIqIw0uZgZ7VaYbVaQ/LiQ4cOxdy5c1FcXIyUlBQAwGeffQa9Xo/BgweH5DWIwoWhTx9krHgFeRNuQf327Si44w5kLF8GlZnhjoiIjk27HmOXn5+P7du3Iz8/H36/H9u3b8f27dthb7hW5kUXXYR+/fph7Nix2LZtG7788kvMnDkTkyZN4hmu1CUZTj4ZGS+/DFV0NOq3bEHB5DshOp1Kt0VERGGiXYc7mTBhAlatWnXE8g0bNmDkyJEA5PA3ZcoUfPXVVzAajbjxxhvx9NNPH/PuVg53QpGo/rffkH/rRIh2O0xnn430pUugMhiUbouIiBTQlqzTIePYtScGO4pUzm3bUDDxNohOJ8znnYv0hQsh6HRKt0VERB0sbMaxI6LWmQYNQvryZRAMBji+/Q4H77sPkt+vdFtERNSJMdgRdWKmIUOQtmABoNWi7uNPUDJnDsJ8IzsREbUjBjuiTi7q3OHo/tRTgEqFmrf+i7L5TzHcERFRixjsiMJAzMV/RcpjjwEAqlasQMWSJQp3REREnRGDHVGYiL3maiQ9+AAAoOKFBah69TWFOyIios6GwY4ojMSPGwfrXdMAAKWPP46ad95VuCMiIupMGOyIwox1yhTEjx8PACh+6CHUfvqZwh0REVFnwWBHFGYEQUC3+++DZcw1gCji4MyZsH/3vdJtERFRJ8BgRxSGBEFAyiOPIPqSiwGvF4V33YX67duVbouIiBTGYEcUpgS1Gt2ffBLm886F5HKh4I7JcB84oHRbRESkIAY7ojAm6HRIe/55GE45BX6bDQW3TYK3tEzptoiISCEMdkRhTmUyIX3ZUuiysuAtKkLB7bfDX1endFtERKQABjuiCKCJi0P6Sy9CnWiFe88eFE6ZCtHtVrotIiLqYAx2RBFCl5aGjOXLoTKb4dy8GUX/ug+S3690W0RE1IEY7IgiiOHkk5G2aCEErRZ1n36K0sfn8bqyRERdCIMdUYQxn302Uuc/CQgCql9/HZUvvqR0S0RE1EEY7IgiUMwllyDpAfm6suXPPstLjxERdREMdkQRKn7cWCRMug0AUPzww7B/843CHRERUXtjsCOKYIn33APLlVcCfj8KZ/wD9Tt2Kt0SERG1IwY7oggmCAJS/v0YzOecA6m+HgV3Toa3qEjptoiIqJ0w2BFFOEGrRff/PA99797wl1egYPKd8NvtSrdFRETtgMGOqAtQR0UhfekSeQDjvXtxcMY/IHm9SrdFREQhxmBH1EVoU1ORvmQpBKMRju+/R8m/53KMOyKiCMNgR9SFGAf0R/dnngYEATVr16LqlRVKt0RERCHEYEfUxUSffz6SHrgfAFD21FOo/fQzhTsiIqJQYbAj6oLixo5F3E03AQCK/vUv1P/6q8IdERFRKDDYEXVBgiAg6cEHEDVyJCS3GwVTpsJTWKh0W0REdIIY7Ii6KEGtRvdnnoa+38nwV1ai4I7J8NtsSrdFREQngMGOqAtTmc1IX7IEmqQkeP78E4UzZnAYFCKiMMZgR9TFaZOSkL5sKVQmE5w//YzSeU8o3RIRER0nBjsigqFvX6Q+NR8QBFS/8Qaq16xRuiUiIjoODHZEBACIvuACJM6YAQAo+fdcODZuUrYhIiJqMwY7IgpIuH0SYi67DPD5cHD6dHgKCpRuiYiI2oDBjogCBEFAyr8fg2HgQPhtNhROmQK/3a50W0REdIwY7IioGZXBgLSFC6Hp1g3ufftRNPOfkPx+pdsiIqJjwGBHREfQJnVD2qKFEPR62L/+GuXPP690S0REdAwY7IioRcaBA5Eydy4AoPLFl2B7/32FOyIioqNhsCOiVlkuuxQJd9wBACh++P9Qv327sg0REVFQDHZEFFTi3dMRdcEFkDweFEy7C97iYqVbIiKiVjDYEVFQgkqF7vOfhL53b/grKlB413SIbrfSbRERUQsY7IjoqFRmM9IWL4Y6NhauHTtQMucRSJKkdFtERHQYBjsiOia6tO7o/uwzgEoF27vvovqNN5RuiYiIDsNgR0THzDxsGLrNnAkAKJ33BJy//KJwR0RE1BSDHRG1SfwtExAzejTg86Hw7hnwlpQo3RIRETVgsCOiNmm87Ji+Tx/4KytROP1unkxBRNRJMNgRUZupTCakLVoItcUC12+/oeTRR3kyBRFRJ8BgR0THRZeWhtTGkynefgc1a9cq3RIRUZfHYEdExy3qnHPQ7Z5/AABK5j4O59atCndERNS1MdgR0QmJnzgR0ZdcDHi9KJx+N7ylpUq3RETUZTHYEdEJEQQBqXPnBq5McXD63RA9HqXbIiLqkhjsiOiEqUwmpC1cAJXFgvpff0Xpv+cq3RIRUZfEYEdEIaHLyED3p58GBAE169ahet06pVsiIupyGOyIKGSizh2OxLvvBgCUPvZv1P/6q8IdERF1LQx2RBRSCXfcjugLL4TUcDKFr6JC6ZaIiLoMBjsiCilBEJAybx50PXvCV1qKwhkzIHm9SrdFRNQlMNgRUcipo8xIW7AAqqgo1P+yBaVPzle6JSKiLoHBjojahT67B1LnPwkAqF69Grb331e4IyKiyMdgR0TtJvr882GdMgUAUPx/s1G/c6fCHRERRTYGOyJqV9ZpUxE1YgQktxsH75oOX3W10i0REUUsBjsialeCSoXUp+ZDm5kBb1ERDt5zDySfT+m2iIgiUrsGu7lz52LYsGEwmUyIjY1tsUYQhCOmpUuXtmdbRNTB1DExSFuwAILJBOdPP6PsueeUbomIKCK1a7DzeDy49tprceeddwatW7FiBYqLiwPT+PHj27MtIlKAoXdvpM79NwCg6uVXYFv/ocIdERFFHk17rvyRRx4BAKxcuTJoXWxsLJKTk9uzFSLqBGIuuQSuXbtQ+eJLKJ41C7qsLBgH9Fe6LSKiiNEpjrGbNm0arFYrzjjjDCxduhSiKCrdEhG1k8QZM2AecR4ktxuF06bBV16udEtERBFD8WD32GOP4a233sIXX3yB66+/Hvfeey8ef/zxVuvdbjdqa2ubTUQUPgS1Gt2ffhq6Hj3gKylB4fS7IXo8SrdFRBQR2hzs5syZ0+IJD02nX3755ZjX99BDD2Ho0KE47bTTcO+99+LRRx/FU0891Wr9vHnzYLFYAlN6enpbfwQiUpg6OhppixdBFR2N+m3bUPLoo5AkSem2iIjCniC18dO0oqICFUe5qHdWVhYMBkPg/sqVKzFjxgzU1NQcdf0//PADhg8fjpKSEiQlJR3xuNvthtvtDtyvra1Feno6bDYbYmJijv0HISLF2b/7DgV3TAZEEUmzZiF+7M1Kt0RE1OnU1tbCYrEcU9Zp88kTVqsVVqv1uJs7mm3btsFgMLQ6PIper4der2+31yeijhN17rnodu+9KHvqKZQ+8QT0J/WEeehQpdsiIgpb7XpWbH5+PqqqqpCfnw+/34/t27cDAE466SRERUXhgw8+QElJCYYOHQqj0YgNGzZg1qxZuP322xneiLqI+FtvgWvPbtT+7wMcnPEPZP33Leh4iAUR0XFp867YtpgwYQJWrVp1xPINGzZg5MiR+OSTT/DAAw9g//79EEUR2dnZuO222zB16lRoNMeWOduyeZKIOifR5ULe2HFw/f479L1OQuaba6COMivdFhFRp9CWrNOuwa4jMNgRRQZvaSlyx1wLX3k5oi64AGkLXoCgUvzEfSIixbUl6/BTk4g6BW1SkhzmtFrYv/wSFQsXKt0SEVHYYbAjok7DeNppSH70UQBAxeIlsH3wgcIdERGFFwY7IupUYv92FeIn3goAKH5wFpxtGBeTiKirY7Ajok6n2733IvrCCyF5vSicdhc8eXlKt0REFBYY7Iio0xFUKqTOfxKGgQPhr6lBwe13wFddrXRbRESdHoMdEXVKKqMR6YsXQZOaAk9eHg7eNZ3XlCUiOgoGOyLqtDSJiUhfuhSqqCg4f/kFJQ8/zGvKEhEFwWBHRJ2aoXdvdH/+eUCthu39/6FiyRKlWyIi6rQY7Iio04safg6SH34YAFDxwgLY1n+ocEdERJ0Tgx0RhYW46/+O+FsbhkF54AE4t2xRuCMios6HwY6Iwka3mfci+sK/yMOgTJ3GYVCIiA7DYEdEYUMeBmU+DAMGwF9Tg/zbb4evslLptoiIOg0GOyIKKyqjEWmLF0HbvTu8efkouP0O+O0OpdsiIuoUGOyIKOxou3VD+ksvQh0XB9fOnTg4/S6OcUdEBAY7IgpT+h49kL58GQSTCY4ff0Lx/fdDEkWl2yIiUhSDHRGFLePAgUhb8AKg1aL2o49ROvdxDmBMRF0agx0RhbWoc85B6rx5AIDq119H5bJlCndERKQcBjsiCnuWyy5F0oMPAgDKn/8PqtetU7gjIiJlMNgRUUSIHzcWCXfcAQAomfMIaj//XOGOiIg6HoMdEUWMxBl3wzLmGkAUUXTvTDg3b1a6JSKiDsVgR0QRQxAEpMyZg6jzz4fk8aBgylS4du9Wui0iog7DYEdEEUXQaND92WdgHDwYYl0d8m+dCPf+/Uq3RUTUIRjsiCjiqAwGpC9ZDEO/fvBXVSHvllvgzslRui0ionbHYEdEEUkdE4P0l1+Cvndv+MsrkD/hFngKCpRui4ioXTHYEVHE0sTFIWPFK9D17AlfaSnyx0+A9+BBpdsiImo3DHZEFNE0CQlyuMvMhLeoCHm33ApvaanSbRERtQsGOyKKeNpu3ZCxaiW0aWnw5ucjf8It8JWXK90WEVHIMdgRUZegTU5GxsqV0KSmwJOTg/xbb4WvqkrptoiIQorBjoi6DF1ad2SuXAlNt25w79uP/Fsnwl9To3RbREQhw2BHRF2KLiMDGStXQm21wr17N/JvmwR/ba3SbRERhQSDHRF1OfrsHshc8QrUcXFw7dghH3NXXa10W0REJ4zBjoi6JH2vXshYuQLq+Hi4du1C/rhx8JaVKd0WEdEJYbAjoi7L0KcPMle/duiYu7Hj4C0uVrotIqLjxmBHRF2aPjsbma+vhrZ7d3jy8pB3083w5Ocr3RYR0XFhsCOiLk+Xno7M1a9Bl5UlD2J8081w//mn0m0REbUZgx0REQBtSgoyX3sV+l694CsvR97YcXDt3q10W0REbcJgR0TUQJOYiIxXV8HQvz/8VVXIGzce9b/+qnRbRETHjMGOiKgJTVwcMlaugHHQIIi1tci/5VY4N29Wui0iomPCYEdEdBh1dDQyXn4JprPPhuh0In/S7aj76iul2yIiOioGOyKiFqhMJqQvXYKokSMhuVwonHYXqtesVbotIqKgGOyIiFqhMhiQtnABLGOuAUQRJXPmoPyFFyBJktKtERG1iMGOiCgIQaNBymOPwTp1KgCgYvESFM96CJLXq3BnRERHYrAjIjoKQRCQeNc0JD/2KKBWw/bOOyiYMhWiw6F0a0REzTDYEREdo7hrr0XawgUQDAY4vvsOeePGw1dRoXRbREQBDHZERG0QPWoUMlethDouDq6dO5F7/Q1w5+Qo3RYREQAGOyKiNjOeeiqy3nwD2vR0eAsLkXfDjajfvl3ptoiIGOyIiI6HLisLWW++AcOAAfDX1CBv/ATY1n+odFtE1MUx2BERHSeN1YrMVSvlse7cbhTNnImy556HJIpKt0ZEXRSDHRHRCVCZzUhbtBAJt00EAFQuW4bCu6bDb+cZs0TU8RjsiIhOkKBWo9vMmUid/yQEnQ72L79E3g03wFNYqHRrRNTFMNgREYWI5YorkPnaq1AnWuHetw+5Y66FY9Mmpdsioi6EwY6IKISMp56KHm+9BUP//vDX1CD/1omoXrtO6baIqItgsCMiCjFtcjIyX1+NmNGjAZ8PJbNno+Sxf/MyZETU7hjsiIjagcpgQOozTyNxxt0AgOrXX0f+xNvgKy9XuDMiimQMdkRE7UQQBFgnT0baooVQmUxwbtqEA1dfzePuiKjdMNgREbWz6AsuQNZ//wt9r5PgL69A/oRbULH8RY53R0Qhx2BHRNQB9Nk9kLV2LSxXXgGIIsqffRaFU6bCX1OjdGtEFEEY7IiIOojKZELKE08g+dFH5PHuvv4aOVdfg/rff1e6NSKKEAx2REQdSBAExF13HbLWvAltRga8RUXIu/EmVL3xBiRJUro9IgpzDHZERAow9OuHHv99C1F/uQCS14vSRx9D0b0zeSkyIjohDHZERApRx8QgbcECdLvvPkCtRu1HHyHn6qtR/+uvSrdGRGGq3YJdbm4uJk6ciB49esBoNKJnz56YPXs2PB5Ps7r8/HxcfvnlMJvNsFqtmD59+hE1RESRShAEJNwyAZmvvQpNSgq8+fnIvfEmlC9aBMnnU7o9Igoz7Rbsdu/eDVEUsWzZMuzcuRPPPfccli5digcffDBQ4/f7cemll8LhcOD777/HmjVr8Pbbb+Pee+9tr7aIiDol0+mnI/v99xBz6aWA34+KBQuRd/NYeAoKlG6NiMKIIHXg0bpPPfUUlixZggMHDgAAPv74Y1x22WUoKChAamoqAGDNmjWYMGECysrKEBMTc9R11tbWwmKxwGazHVM9EVFnZ/vgA5Q88ihEux0qkwlJDz0Ey9+ugiAISrdGRApoS9bp0GPsbDYb4uPjA/d/+uknDBgwIBDqAOCvf/0r3G43tmzZ0uI63G43amtrm01ERJHEcvnl6PHeezAOGQzR6UTxgw/i4Ix/cMw7IjqqDgt2f/75JxYsWIDJkycHlpWUlCApKalZXVxcHHQ6HUpKSlpcz7x582CxWAJTenp6u/ZNRKQEXVp3ZK5ahcQZMwCNBnWffooDV14Fx08/Kd0aEXVibQ52c+bMgSAIQadffvml2XOKiopw8cUX49prr8Vtt93W7LGWdi1IktTqLocHHngANpstMBXw+BMiilCCWg3r5DuQ9eab0GVlwVdaivxbbkXJY/+G6OCwKER0JE1bnzBt2jRcf/31QWuysrICt4uKijBq1CgMHToUy5cvb1aXnJyMjRs3NltWXV0Nr9d7xJa8Rnq9Hnq9vq1tExGFLePAAejxztsofXI+atauRfXrr8O+YQOSH3sUUeeco3R7RNSJtOvJEwcPHsSoUaMwePBgrF69Gmq1utnjjSdPFBYWIiUlBQCwdu1ajB8/nidPEBG1wPHjjyh++P/gPXgQAGC55mok3Xcf1Pz8I4pYbck67RbsioqKMGLECGRkZODVV19tFuqSk5MByMOdnHbaaUhKSsJTTz2FqqoqTJgwAVdddRUWLFhwTK/DYEdEXY3ocKDsuedR/frrgCRBk5iI5EfmIPr885VujYjaQacIditXrsQtt9zS4mNNXzI/Px9TpkzBV199BaPRiBtvvBFPP/30Me9uZbAjoq7KuXUrih+cBU9uLgAgZvRoJD00C5omow8QUfjrFMGuozDYEVFXJrpcqFi0CJWvrAD8fqjj4pA0axZiLh3Nce+IIkSnHceOiIhCS2UwoNu99yJr7Vro+/SBv7oaRTNnomDibXAfyFG6PSLqYAx2REQRwDigP3q8tQ7W6XdB0Ong+PFHHLjySpQ99zzE+nql2yOiDsJgR0QUIQSdDolTpiB7/QcwjzgP8HpRuWwZDlx6Geq+/BJhfuQNER0DBjsiogijy8hA+tKlSFu4AJrUFHiLilA4dRoKJ98JDwd1J4poDHZERBFIEARE/+Uv6Ll+PRJuvx3QamH/5hscuOxylC9aBNHtVrpFImoHDHZERBFMZTKh2z3/QPb778E09GxIbjcqFizEgUsvQ+0nn3L3LFGEYbAjIuoC9NnZyHjlFXR/9hlounWDt7AQB2fMQN5NN6P+11+Vbo+IQoTBjoioixAEATGjR6Pnxx/BOnUqBKMR9Vu3Ivfv1+PgvTMDlykjovDFYEdE1MWozGYk3jUNPT/5GJa//Q0QBNR++CH+vGQ0yp55Fn67XekWieg4MdgREXVR2qQkpM57HD3e/i9MZ50FyeNB5Ysv4s+L/orqNWsg+XxKt0hEbcRLihERESRJgn3DBpTNfypw7Vldjx5IvGsaoi++GIKK2wGIlMJrxRIR0XGRvF5Ur1mLikWL4K+pAQDo+/ZF4vTpiBo1ktefJVIAgx0REZ0Qv92OqpWrULViBUSHAwBgOPUUdJsxA+ahQxXujqhrYbAjIqKQ8FVXo+qVV1D12mpILhcAwHTWWUi8+26YTh+kcHdEXQODHRERhZSvvBwVy5ajZu1aSF4vAMA84jwkTpsG48CBCndHFNkY7IiIqF14i4pQsWQJat55F/D7AQDmYcOQMPkOmM44g8fgEbUDBjsiImpXntxcVCxdBtsHHwQCnvH002GdfAfM557LgEcUQgx2RETUITyFB1H58kuwvf0OJI8HAGDo1w8Jd9yB6Av/wmFSiEKAwY6IiDqUt6wMVStWonrtWkhOJwBA17MnrLdPQszo0RC0WoU7JApfDHZERKQIX3U1ql97DVWrX4dYWwsA0CQlIe7mmxB33XVQWywKd0gUfhjsiIhIUX67HdVvvImq116Fv7wCACAYjYi9+mrEjxsLXWamwh0ShQ8GOyIi6hREjwe1H36EqpUr4d6zR14oCIg6/3wkTBgP45AhPNGC6CgY7IiIqFORJAnOjRtRtWIl7N98E1hu6N8f8ePHIfrii6HS6RTskKjzYrAjIqJOy33gAKpWvQrbe+9BcrsBAOr4eMReczVi//536NLSFO6QqHNhsCMiok7PV12NmrVrUb1mLXwlJfJCQYD5vHMRd/31iDrvPAhqtbJNEnUCDHZERBQ2JJ8P9q+/RvWba+D44YfAcm1qKmL//nfEjrkGmoQEBTskUhaDHRERhSVPbi6q166D7Z134LfZ5IVaLWIu/AssV18D89CzuRWPuhwGOyIiCmuiy4XaTz5BzZtrUP/rr4HlmuRkWK66ErF/+xuHTKEug8GOiIgihuuPP1Dz9juwffABxMateABMQ4bAcvXViPnrRVCZzQp2SNS+GOyIiCjiiB4P7F99hZq334Hj+++Bhq8vlcmE6EsuhuXKK2EaMoTXp6WIw2BHREQRzVtSAtt776Pm3XfgzcsPLNckJyPm0tGwXHYZ9H37cvBjiggMdkRE1CVIkoT6LVtQ8+67qPvsc4h1dYHHdD17wnLZpYi57DLo0tMV7JLoxDDYERFRlyO63bB/+y1q138I+4YNkDyewGPGU09FzGWXIfqii6BN6qZgl0Rtx2BHRERdmr+uDnWff4Ha9evh+PlnQBTlBwQBxkGDEH3RhYi58EJou3dXtlGiY8BgR0RE1MBXXo7ajz9B7YcfNhs6BQAMAwYg+qKLEHPRhdBlZSnTINFRMNgRERG1wFtSgrrPv0DdZ5/BuWXLoS15APR9+iD6wgsRff4o6E8+mSdeUKfBYEdERHQUvspK1H3xJeo+/RSOjRsBvz/wmCY5GVEjRyD6/PNhOussqPR6BTulro7BjoiIqA38NTWo+2oD6r78Eo4ff4RUXx94TDCZYB42FNGjRiFqxAhorFYFO6WuiMGOiIjoOIkuF5wbN6JuwwbYN3wNX2npoQcFAYaBAxE1/ByYhw+H8ZRTIGg0yjVLXQKDHRERUQhIkgT3H38EQp5rx45mj6uio2E++2yYzx2OqOHDoU1NVahTimQMdkRERO3AW1oGx/ffw/HD93D88CP8Ta5dCwC67GyYh58D87BhMA0ZAnVUlEKdUiRhsCMiImpnkt8P186dsH//PRzffY/6335rdgIG1GoYBwyA6ayzYD77LBgHDYLKaFSuYQpbDHZEREQdzF9bC8dPP8tb9DZuhDc/v9njglYL46mnwnT22TCfdSYMp54KlU6nULcUThjsiIiIFOYtKoJj4yY4f/4Zjo0b4Sspafa4oNPBcMpAmE4fDNOQwTAOGgR1dLRC3VJnxmBHRETUiUiSBG9+Phw/b4Rz489wbNwEf2Vl8yJBgL5PH5gGD4Zp8OkwDh7C69oSAAY7IiKiTk2SJHhyc1G/dSucv2yBc8uWI3bdAoAmNQXGU08NTIZ+/ThYchfEYEdERBRmvGVlctDbshXOLb/AvXtPs0ueAQC0WhhOPrlJ2DsF2rQ0Xv4swjHYERERhTm/3QHXjh2o//XXwHTE7lsAKosFxv79YejfH4YBA2Do3x/a7qkMexGEwY6IiCjCSJIE78GDqN/eEPR++xXuXX9A8nqPqFXHxjYJev1g6NtX3rKnUinQOZ0oBjsiIqIuQPJ44Nq3D64dO+HasQOunTvh2rsX8PmOqFWZzdD37QtDnz7Qn9wXhr59oe/VCyqDQYHOqS0Y7IiIiLoo0e2Ge+9euHbuRP2OHXDv+gPuffta3LIHlQq67B4w9O4Dfa+ToO/VC/peveSte2p1xzdPLWKwIyIiogDJ64U7Jwfu3bvh2r0H7t1/wPXHbvirq1usFwwG6Hv2bAh6DYGvZ09oUlK4O1cBDHZEREQUlCRJ8JWVw71b3qInT/vh/vNPSG53i88RjEboemRB3yMbup7Z0GdnQ5edDV1WFq+i0Y4Y7IiIiOi4SH4/vAUFcDWEPc/+/XLoy80DWtqdCwAqFbTpadBlZUGXmdkwZUGXlQVtSjJ3654gBjsiIiIKKcnng6egAJ6cHLj//BOeAzlwH/gTnj8PQLTbW32eoNVCm5FxKPBlpEObli7PU1IgcEvfUTHYERERUYeQJAm+8nJ4DuTAk5cnT7m58OTlwZuf3/JJG41UKmhTUqBNT4cuPb1hngZt9+7QpqZCnZDA8fjAYEdERESdgOT3w1tcDE9uHjx5DWGvoBDewgJ4CgohuVxBny8YDNCmpgaCnrZ7d2i7p0KbkgptSjI0iYkQNJoO+mmUw2BHREREnVrjlj5vQQE8BQXwFhTCU5APb+FBeIuK4CstBY4WUdRqaLp1gzY5WQ56KSnQJqfIt5OSoUnqBk1CQtgf48dgR0RERGFN8njgLSmBt6gI3oMHG6YieA4WwldcAm9paYsDMR9BrYYmMRGapG7QdkuCJilJvp2UJC9vmFQxMZ12t29bsk7kb78kIiKisCPodNBlZECXkdHi45LfD19FJXwlxfAWl8BbXNz8dmkpfBUVgN8PX0kJfCUlCLbjV9DrobFam4U9TaIVaqsVmgQrNNYEaBISoLZaodLr2+eHDoF2C3a5ubl47LHH8NVXX6GkpASpqam4+eabMWvWLOianAHTUjpesmQJJk+e3F6tERERUZgT1Gpok7pBm9QNxlNPbbFG8vngq6yEr7QU3tJS+ErL5MBXVgpvaRl85eXwlZdDrK2F5HYHtgwejSoqKhDyNAkJSH74IWgSE0P9Ix6Xdgt2u3fvhiiKWLZsGU466STs2LEDkyZNgsPhwNNPP92sdsWKFbj44osD9y0WS3u1RURERF2EoNFAm5QEbVISjEHqRJcLvooK+MrKA2HPV14OX0U5/BWVcjisrIS/ogKS1wvRbofHbgfy8gAAyY/M6ZCf51i0W7C7+OKLm4W17Oxs7NmzB0uWLDki2MXGxiI5Obm9WiEiIiJqlcpggC4tDbq0tKB1kiRBrKuDr6IS/soKOQxWVELdiTZIdegxdjabDfHx8UcsnzZtGm677Tb06NEDEydOxO233w5VK9eic7vdcDe51EltbW279UtERETUSBAEqGNioI6JAbJ7KN1Oizos2P35559YsGABnnnmmWbLH3vsMVxwwQUwGo348ssvce+996KiogIPPfRQi+uZN28eHnnkkY5omYiIiCistHm4kzlz5hw1WG3evBlDhgwJ3C8qKsKIESMwYsQIvPTSS0Gf+8wzz+DRRx+FzWZr8fGWttilp6dzuBMiIiKKSO06jl1FRQUqKiqC1mRlZcFgMACQQ92oUaNw1llnYeXKla3uYm30ww8/YPjw4SgpKUFSUtJR++E4dkRERBTJ2nUcO6vVCqvVeky1Bw8exKhRozB48GCsWLHiqKEOALZt2waDwYDY2Ni2tkZERETUpbXbMXZFRUUYOXIkMjIy8PTTT6O8vDzwWOMZsB988AFKSkowdOhQGI1GbNiwAbNmzcLtt98OfSce/I+IiIioM2q3YPfZZ59h//792L9/P9IOO324ce+vVqvF4sWLcc8990AURWRnZ+PRRx/F1KlT26utE+JwOFp9TK1WB3Y/H61WpVLBaDQeV63T6URre88FQYDJZDqu2vr6eoii2GofZrP5uGpdLhf8fn9Iak0mU2BAa7fbDV+QS8m0pdZoNAa2Jns8Hni93pDUGgwGqBuuT9iWWq/XC4/H02qtXq+HpuGi122p9fl8zY5PPZxOp4NWq21zrd/vhyvIhby1Wm1gUPK21IqiiPr6+pDUajSawB+LkiTB6XSGpLYt/+75GdFyLT8j+BkR7p8Rne4yZFKYs9lsEgDJZrO1+2sBaHUaPXp0s1qTydRq7YgRI5rVWq3WVmuHDBnSrDYzM7PV2n79+jWr7devX6u1mZmZzWqHDBnSaq3Vam1WO2LEiFZrTSZTs9rRo0cHfd+aGjNmTNBau90eqB0/fnzQ2rKyskDtlClTgtbm5OQEamfOnBm0dseOHYHa2bNnB63dtGlToHb+/PlBazds2BCoXbhwYdDa9evXB2pXrFgRtHbdunWB2nXr1gWtXbFiRaB2/fr1QWsXLlwYqN2wYUPQ2vnz5wdqN23aFLR29uzZgdodO3YErZ05c2agNicnJ2jtlClTArVlZWVBa8ePHx+otdvtQWvHjBnT7Hc4WC0/I+SJnxGHJn5GyFO4f0Z0hLZknaMf9EZEREREYaHNZ8V2Nh15Vix3s7S9lrtZuJsl3HezcFesjJ8R/IzgZ0TLtR2xK7ZdhzvpbDjcCREREUWytmQd7oolIiIiihAMdkREREQRgsGOiIiIKEIw2BERERFFCAY7IiIiogjBYEdEREQUIRjsiIiIiCIEgx0RERFRhGCwIyIiIooQDHZEREREEYLBjoiIiChCMNgRERERRQgGOyIiIqIIwWBHREREFCE0SjdwoiRJAgDU1tYq3AkRERFR6DVmnMbME0zYB7u6ujoAQHp6usKdEBEREbWfuro6WCyWoDWCdCzxrxMTRRFFRUWIjo6GIAjt9jq1tbVIT09HQUEBYmJi2u11uiq+v+2L72/74vvbvvj+ti++v+0rFO+vJEmoq6tDamoqVKrgR9GF/RY7lUqFtLS0Dnu9mJgY/uK3I76/7Yvvb/vi+9u++P62L76/7etE39+jbalrxJMniIiIiCIEgx0RERFRhGCwO0Z6vR6zZ8+GXq9XupWIxPe3ffH9bV98f9sX39/2xfe3fXX0+xv2J08QERERkYxb7IiIiIgiBIMdERERUYRgsCMiIiKKEAx2RERERBGCwe4YLF68GD169IDBYMDgwYPx3XffKd1SxJg3bx7OOOMMREdHo1u3brjqqquwZ88epduKSPPmzYMgCJgxY4bSrUSUgwcP4uabb0ZCQgJMJhNOO+00bNmyRem2wp7P58NDDz2EHj16wGg0Ijs7G48++ihEUVS6tbD17bff4vLLL0dqaioEQcB7773X7HFJkjBnzhykpqbCaDRi5MiR2LlzpzLNhqFg76/X68V9992HgQMHwmw2IzU1FePGjUNRUVHI+2CwO4q1a9dixowZmDVrFrZt24Zzzz0Xl1xyCfLz85VuLSJ88803mDp1Kn7++Wd8/vnn8Pl8uOiii+BwOJRuLaJs3rwZy5cvxymnnKJ0KxGluroa55xzDrRaLT7++GPs2rULzzzzDGJjY5VuLew9+eSTWLp0KRYuXIg//vgD8+fPx1NPPYUFCxYo3VrYcjgcOPXUU7Fw4cIWH58/fz6effZZLFy4EJs3b0ZycjIuvPDCwDXZKbhg76/T6cTWrVvx8MMPY+vWrXjnnXewd+9eXHHFFaFvRKKgzjzzTGny5MnNlvXt21e6//77FeoospWVlUkApG+++UbpViJGXV2d1KtXL+nzzz+XRowYId19991KtxQx7rvvPmn48OFKtxGRLr30UunWW29ttuzqq6+Wbr75ZoU6iiwApHfffTdwXxRFKTk5WXriiScCy1wul2SxWKSlS5cq0GF4O/z9bcmmTZskAFJeXl5IX5tb7ILweDzYsmULLrroombLL7roIvz4448KdRXZbDYbACA+Pl7hTiLH1KlTcemll+Ivf/mL0q1EnP/9738YMmQIrr32WnTr1g2DBg3Ciy++qHRbEWH48OH48ssvsXfvXgDAr7/+iu+//x6jR49WuLPIlJOTg5KSkmbfd3q9HiNGjOD3XTux2WwQBCHkW/g1IV1bhKmoqIDf70dSUlKz5UlJSSgpKVGoq8glSRLuueceDB8+HAMGDFC6nYiwZs0abN26FZs3b1a6lYh04MABLFmyBPfccw8efPBBbNq0CdOnT4der8e4ceOUbi+s3XfffbDZbOjbty/UajX8fj/mzp2LG264QenWIlLjd1pL33d5eXlKtBTRXC4X7r//ftx4442IiYkJ6boZ7I6BIAjN7kuSdMQyOnHTpk3Db7/9hu+//17pViJCQUEB7r77bnz22WcwGAxKtxORRFHEkCFD8PjjjwMABg0ahJ07d2LJkiUMdido7dq1WL16Nd544w30798f27dvx4wZM5Camorx48cr3V7E4vdd+/N6vbj++ushiiIWL14c8vUz2AVhtVqhVquP2DpXVlZ2xF81dGLuuusu/O9//8O3336LtLQ0pduJCFu2bEFZWRkGDx4cWOb3+/Htt99i4cKFcLvdUKvVCnYY/lJSUtCvX79my04++WS8/fbbCnUUOf75z3/i/vvvx/XXXw8AGDhwIPLy8jBv3jwGu3aQnJwMQN5yl5KSEljO77vQ8nq9uO6665CTk4Ovvvoq5FvrAJ4VG5ROp8PgwYPx+eefN1v++eefY9iwYQp1FVkkScK0adPwzjvv4KuvvkKPHj2UbiliXHDBBfj999+xffv2wDRkyBDcdNNN2L59O0NdCJxzzjlHDM+zd+9eZGZmKtRR5HA6nVCpmn9FqdVqDnfSTnr06IHk5ORm33cejwfffPMNv+9CpDHU7du3D1988QUSEhLa5XW4xe4o7rnnHowdOxZDhgzB0KFDsXz5cuTn52Py5MlKtxYRpk6dijfeeAPvv/8+oqOjA1tHLRYLjEajwt2Ft+jo6COOVTSbzUhISOAxjCHyj3/8A8OGDcPjjz+O6667Dps2bcLy5cuxfPlypVsLe5dffjnmzp2LjIwM9O/fH9u2bcOzzz6LW2+9VenWwpbdbsf+/fsD93NycrB9+3bEx8cjIyMDM2bMwOOPP45evXqhV69eePzxx2EymXDjjTcq2HX4CPb+pqamYsyYMdi6dSvWr18Pv98f+L6Lj4+HTqcLXSMhPcc2Qi1atEjKzMyUdDqddPrpp3MojhAC0OK0YsUKpVuLSBzuJPQ++OADacCAAZJer5f69u0rLV++XOmWIkJtba109913SxkZGZLBYJCys7OlWbNmSW63W+nWwtaGDRta/LwdP368JEnykCezZ8+WkpOTJb1eL5133nnS77//rmzTYSTY+5uTk9Pq992GDRtC2ocgSZIUuphIRERERErhMXZEREREEYLBjoiIiChCMNgRERERRQgGOyIiIqIIwWBHREREFCEY7IiIiIgiBIMdERERUYRgsCMiIiKKEAx2RERERBGCwY6IiIgoQjDYEREREUUIBjsiIiKiCPH//puRd4SGe0kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the results\n",
    "# plt.subplot(2, 1, 1)\n",
    "for i, y in enumerate(C @ xout):\n",
    "    plt.plot(tout, y)\n",
    "    plt.plot(tout, yd[i] * np.ones(tout.shape), 'k--')\n",
    "plt.title('outputs')\n",
    "\n",
    "# plt.subplot(2, 1, 2)\n",
    "# plt.plot(t, u);\n",
    "# plot(np.range(Nsim), us*ones(1, Nsim), 'k--')\n",
    "# plt.title('inputs')\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "# Print the final error\n",
    "xd - xout[:,-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
